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ABSTRACT 

We present a new method for the analysis of images, a fundamental task in observa- 
tional astronomy. It is based on the linear decomposition of each object in the image 
into a series of localised basis functions of different shapes, which we call 'Shapelets'. 
A particularly useful set of complete and orthonormal shapelets is that consisting of 
weighted Hermite polynomials, which correspond to perturbations around a circular 
gaussian. They are also the eigenstates of the 2-dimensional Quantum Harmonic Os- 
cillator, and thus allow us to use the powerful formalism developed for this problem. 
Among their remarkable properties, they are invariant under Fourier transforms (up 
to a rescaling), leading to an analytic form for convolutions. The generator of linear 
transformations such as translations, rotations, shears and dilatations can be written 
as simple combinations of raising and lowering operators. We derive analytic expres- 
sions for practical quantities, such as the centroid (astrometry) , flux (photometry) 
and radius of the object, in terms of its shapelet coefficients. We also construct polar 
basis functions which are eigenstates of the angular momentum operator, and thus 
have simple properties under rotations. As an example, we apply the method to Hub- 
ble Space Telescope images, and show that the small number of shapelet coefficients 
required to represent galaxy images lead to compression factors of about 40 to 90. 
We discuss applications of shapelets for the archival of large photometric surveys, for 
weak and strong gravitational lensing and for image deprojection. 

Key words: methods: data analysis, analytical; techniques: image processing, sur- 
veys, gravitational lensing 



1 INTRODUCTION 

Image analysis is a fundamental task in observational as- 
tronomy. For instance, new techniques, such as weak grav- 
itational lensing (see reviews by Mellier 1999; Bartelmann 
& Schneider 2000), microlensing (Mao 1999) and the search 
for supernovae (Riess 2000; Perlmutter et al. 1998), have 
great scientific promise, but require high precision image 
processing and analysis. As a result, a number of sophisti- 
cated data analysis packages (eg. FOCAS in IRAF, Jarvis 
& Tyson 1986; SExtractor, Bertin & Arnouts 1996, etc) 
and techniques (eg. wavelet analysis, see review by Stark, 
Murtagh & Bijaoui 1998; image subtraction, Alard & Lup- 
ton 1998; shear measurement, Kaiser, Squires & Broadhurst 
1995, Kaiser 2000, and Kuijken 2000) have been developed. 

In this paper, we present a new method for image anal- 
ysis. It is based on the linear decomposition of each ob- 
ject into a series of localised basis functions with different 
shapes, which we call 'Shapelets'. As a basis set we choose 
weighted hermite polynomials. They correspond to pertur- 
bations about a circular gaussian, and, in their asymptotic 
form, to the Edgeworth expansion in several dimensions. 
They are also the eigenstates of the 2-dimensional Quan- 



tum Harmonic Oscillator (QHO), and thus allow us to use 
the powerful formalism developed for this problem. They 
have remarkable properties. In particular, they are (up to 
a rescaling) invariant under Fourier transforms and thus 
yield a simple analytical form for convolutions. We derive 
a number of practical tools which can use used to compute 
the characteristics of the object (centroid, flux, radius, etc) 
from its shapelet coefficients. Our method differs from the 
wavelet transform which decomposes the image into a sum 
of basis functions of different scales but with a set shape. 
In our method, the image is decomposed into a collection 
of compact disjoint objects of arbitrary shapes and is thus 
particularly adapted to astronomical images. 

As an example, we show how images of galaxies ob- 
served with the Hubble Space Telescope (HST) can be rep- 
resented and strongly compressed using shapelets. We also 
discuss several applications of shapelets, such as archival of 
large photometric catalogues, gravitational lensing and im- 
age de-projection. A precise method to measure the shear 
induced by weak lensing on galaxy images is presented in 
an adjoining paper (Refregier & Bacon 2001, Paper II). The 
application of Shapelets to interferometric images will be 
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Figure 1. First few 1-dimensional basis functions <j> n (x). 



presented in Chang & Refregier (2001). The analytical re- 
sults derived in this paper may also be useful for any ap- 
plication using the Edgeworth expansion, such as, for in- 
stance, the study of the growth of cosmological perturba- 
tions (Juskiewicz et al. 1995 and reference therein). 

This paper is organised as follows. In we describe the 
main properties of 1-dimensional shapelets and discuss their 
connection to the QHO. In we show how 2-dimensional 
shapelets can be formed and derive a number of practical 
analytical results. In we discuss how the shapelet states 
behave under convolutions. In we derive polar shapelets 
from the cartesian basis functions and describe some of their 
properties. In §[], we discuss several direct applications of 
shapelets. Our conclusions are presented in §hj 



2 ONE-DIMENSIONAL SHAPELETS 
2.1 Definitions 

We first consider the description of a localised object in 1- 
dimension. For this purpose, we first define the dimension- 
less basis functions 



<t>n{x) 



2 n 7r2n! 



H n {x)e 2 



(1) 



where n is a non-negative integer and H n (x) is a hermite 
polynomial of order n. These functions are orthonormal in 
the sense that 



dx <)> n (x)<j> m {x) = 5„ 



(2) 



where S mn is the Kronecker delta symbol. The first few func- 
tions are plotted on figure |l[ These functions, which we call 
'Shapelets', can be thought of as shape perturbations around 
the gaussian cj>o(x), 

To describe an object in practice, we use the dimen- 
sional basis functions 



where (3 is a characteristic scale, which is typically chosen 
to be close to the size of the object. These functions are also 
orthonormal, i.e. 



dx B n (x; f3)Bm(x; (3) = <5„„ 



(4) 



This infinite set of functions forms a complete basis for 
smooth and integrable functions. Thus, a (sufficiently well 
behaved) object profile f(x) can be expanded as 



From the orthonormality condition (Eq. 
coefficients are given by 



(•>) 



the shapelet 



dx f(x)B n (x;f3). 



(6) 



In practice, the series of Equation Q will converge quickly 
if the object f(x) is sufficiently localised, and if (3 and the 
origin x — are not too different from the size and loca- 
tion of the object. This series representation is referred to 
as the Gram-Charlier series, or, in its asymptotic form, as 
the Edgeworth expansion (see eg. Juiszkiewicz 1995 and ref- 
erence therein). 

2.2 Fourier Transform 

These basis functions have a number of useful properties. 
Let us first consider their Fourier transform, which, for an 
arbitrary function f(x), is defined as 



/(*) = (2tt)" 
f{x) = (2tt)- 



dxf(x)e l 



dkf(k)e~ 



(7) 



With these conventions, the Fourier transform of the dimen- 
sionless basis function (j> n {£) is 



0„(k) = t n (f>n(K), 



(8) 



Thus, up to a phase factor, the dimensionless basis func- 
tions are invariant under Fourier transforms. This very use- 
ful property can be understood in physical terms fr om the 
analogy with the quantum harmonic oscillator (see §2.3). 



The Fourier transform of the dimensional basis function 
B n {x; (3) is given by 



B n (k;f3)=i n B n (k;f3- 1 )- 



(9) 



(3) 



Thus, the Fourier transform acts on the basis functions with 
an unsurprising change of scale (3 — * f3~ L . 



2.3 Analogy with the Quantum Harmonic 
Oscillator 

As we now discuss, the above basis functions are the eigen- 
states of the Quantum Harmonic Oscillator (QHO), which 
allows us to exploit the readily available formalism devel- 
oped for this problem. Let us consider a QHO with mass m 
and natu ral frequency u. If distances are measured in units 
of U and energies in units of htv, the Hamiltonian for 
this system is 
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H=\[x 2 + p 2 } 



(10) 



where x and p are the position and momentum operators 
respectively. In the ^-representation, they are given by 



„ 1 d 
i ox 



(11) 



and commute as [x,p] = i. As is well known, the basis func- 
tions 4> n {x) are the eigenfunctions of the Hamiltonian, with 



H(p n = ( n + - 



(12) 



Clearly, H is symmetric under a permutation of x and p (see 
Eq. |lfj|), thus explaining the invariance of <j> n under Fourier 
transforms (Eq. pi). 

Of particular practical interest are the lowering and 
raising operators, which are defined as 



1 

71 



(x + ip) . 



V2 



ip) ■ 



(13) 



where ' is the Hermitian conjugate. They commute as 
[a, a*] = 1, and act on the basis functions as 

a<f) n = Vn<j> n -i, a '</> n = Vn + l<f> n +i. (14) 
The Hamiltonian can then be rewritten as H — iV+|, where 



the number operator N = a'a has the property that 
Ncj> n = n<t> n . 



(15) 



When convenient, we will use the bra-ket notation of quan- 
tum mechanics. For instance, the n th state is written as \n) 
and has an ^-representation given by (x\n) = <f> n (x). 

The dimensional basis functions are the eigenfunctions 
of the Hamiltonian 



Hp = \ [P~ 2 x 2 + f3 2 p 2 ] 



(16) 



The eigenstates are labeled as \n; (3) and obey Hp\n;j3) = 
(n + |) [n; /?). The dimensional basis functions are then 
given by B n (x;f3) = (x\n;(3). 



2.4 Further Properties 

The Hermite basis functions have a number of further con- 
venient properties which we will need later and summarise 
here. 

We first notice, by inspecting Figure |l|, that the ba- 
sis functions B n (x,f3) acquire both a larger extent and 
smaller scale oscillations when the order n is increased, keep- 
ing f3 constant. This can be described more precisely by 
considering the characteristic radius S max (/3, n) of a basis 
function, defined by #max(/3,ri) = (n; /3\x 2 \n; (3). As is well 
known from Quantum Mechanics and can easily derived us- 
ing Equation (|l3|), this rms radius is given by 9 mB x(l3,ri) = 

/3(n+i) 5 . Similarly, the characteristic size 9 m i n (/3,n) of 
the small scale (oscillatory) features in a basis function of 
order n is defined as the rms inverse radius in Fourier space, 
i.e. by #~f a (A n) = (n; (3\p 2 \n\ (3) . As can again be verified 
using raising and lower operators, the radius is given by 

9min(P, n) = (3 (n + ^) 5 . Thus a decomposition which in- 
cludes modes from n = to n max can represent features 
with scales ranging between the two limits # m i n (/3, JJ ma x) an d 



#max(/3, ^max). In §3.2 below, we show how these scales can 
be used to fine a good choice of (3 for an object. . 

Another important property relates to the rescaling of 
a shapelet function. Let us for instance consider a function 
f( x ) — X/ fnB n (x; (3), which has been decomposed into 
shapelets of scale /3. It can be sometimes convenient to ex- 
press it in terms of basis functions with a different scale j3' , 
as f(x) = f' n B n (x; f3'). The relation between the coeffi- 
cients /„ and f' n is derived in Appendix [a] and involves the 
overlap m atrix {n; f3\n' , f3') , whose analytic form is given by 
Equation (A3). 

Finally, we note that the basis functions obey the inte- 
gral property 



< l\n;f3 >e 



dx B n (x;{3) 



2 1 - 



n 
n/2 



,(17) 



for n even (the integral vanishes otherwise) , where the paren- 
thesis denotes the binomial coefficient and a convenient 
shorthand notation was used on the left-hand side. This can 
be derived using the generating function of Hermite polyno- 
mials (eg. Arfken 1985). 



3 TWO-DIMENSIONAL CARTESIAN 
SHAPELETS 

In this section, we construct 2-dimensional shapelets by tak- 
ing products of the 1-dimensional shapelets described above. 
We then study the properties of the resulting 'Cartesian' ba- 
sis functions, and derive a number of analytical results which 
are useful in practice. 

3.1 Definitions 

Basis functions for 2-dimensional objects can be constructed 
by taking the tensor product of two 1-dimensional basis 
functions. We thus define the dimensionless functions 



n (x) = 4> ni (xi)(f) n2 (x2), 



(18) 



where x = (x\,X2) and n = (711,712)- Dimensional basis 
functions are defined as. 



(19) 



These 2-dimensional shapelets are again orthonormal, in the 
sense that 



d x B n (x; /3)B m (x; (3) = 5 nimi 5„ 2 



(20) 



The functions </<n(x) are eigenstates of the 2-dimensional 
QHO whose Hamiltonian is 



H 



r-2 ,-2.-2. -21 
\X 1 + X 2 +Pi +p 2 \ 



(21) 



where Xi and pi are the position and momentum operators 
for each dimension. 

The first few 2-dimensional shapelets are shown on Fig- 
ure^ Again, they can be thought of as perturbations around 
the 2-dimensional Gaussian 0oo- These basis functions form 
a complete orthonormal basis for smooth, integrable func- 
tions of two variables. A (well behaved) 2-dimensional func- 
tion /(x), such as the image of an object, can thus be de- 
composed as 
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Figure 2. First few 2-dimensional Cartesian basis functions 
rj> ni ,n 2 ■ The dark and light regions correspond to positive and 
negative values, respectively. 



/(x)= J2 fnBntefl, 



0.0000 0.0005 0.0010 0.0015 0.0020 



original 



ng5 




6 



n^lO 



n^20 





(22) 



Figure 3. Decomposition of a galaxy image found in the HDF. 
The original 60 X 60 pixel HST image (upper left-hand panel) 
can be compared with the reconstructed images with different 
maximum order n = n\ + n 2 . The shapelet scale is chosen to be 
/3 = 4 pixels. The lower right-hand panel (n < 20) is virtually 
indistinguishable from the initial image. 



where the shapelet coefficients are given by 



fn 



d 2 x/(x)S n (x;/3) 



(23) 



Figure ^| show how an image observed with HST can be de- 
composed and reconstructed using shapelets. The resulting 
distribution of the coefficients is shown on Figure ^. More 
examples can be found on Figure |H| These examples and 
associated applications will be discussed in detail in §|(| 

Of practical interest, is the choice of an appropriate 
shapelet scale /3 and maximum order n max for the faithful 
and efficient decomposition of a given image. Using argu- 
ments similar to those of §2.4, it is easy to show that a 



decomposition in 2-dimensions which include shapelets of 
scale (5 with order ranging from n\ + n 2 = to n max can 
only describe features with scales between the two limits 



P (n max + 1)" 



/3(n max + I)'- 



(24) 



Thus, if the function has features with scales ranging from 
# max (eg. the size of the object or that of the image) and 
dmin (eg. the pixel size, or the size of a smoothing kernel), a 
good choice of /3 and n max will be 



f3 S=S (# m in#n 



(25) 



In practice, this provides a good f irst guess, which can be 
refined using a few iterations (see §3.2). 



3.2 Photometry and Astrometry 

The most basic quantities to measure for an object image 
are its total flux (photometry), centroid position (astrom- 
etry) and size. Let us first decompose the intensity /(x) 
of the object into shapelet coefficients / n = (n;/3|/) as in 
Equation Q. 

Using the integral property of Equation (0), it is then 
easy to show that the total flux F = J d 2 xf(x) of the object 
is 



n 2 
n 2 /2 



fn 1 n 2 5 (26) 



where the sum is over even values of m and m. 

Using Equations ( |t^ ) and (E3|), one can also show that 
the centroid of the object x{ = jd 2 xxif(x.)/F is given by 



"2~ 



ni + 1 
(ni + l)/2 



n 2 /2 



(2-ni-n 2 ) 



(27) 



and similarly for x^. 



Similarly, the rms radius r/ defined by rj = 
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■10" 



2 4 6 



10 12 



Shapelet Coefficients 




N Dl = 61x61 



N Dl = 51x51 




N„, =29x29 




N oor =60, Compress = 62 



N cor =55, Compress=47 



N oor =10, Compress = 84 




Figure 4. Shapelet coefficients for the image decomposition of 
the previous figure. Since the coefficient array is sparse, the images 
can be reconstructed from the few first largest coefficients. 



J d 2 xx 2 /(x)/ F is given by 

even 



Figure 5. Reconstruction and compression of three HST galaxy 
images using shapelets. The left-hand column shows the orginal 
images extracted from the HDF and list iVpj x their size in pix- 
els. The right-hand column shows their reconstructed image from 
the N co { largest coefficients (in absolute value) of their shapelet 
decomposition. Because the coefficient matrix is typically sparse, 
a large compression factor Af p ; x /7V co f is achieved. The shapelet 
scale was chosen to be (3 = 4 pixels in all 3 cases. 



1 m/2 ) ( n 2 /2 ' /ni " 2 ' 



(28) 



These expressions can be used, by iteration, to find the op- 
timal centre and scale of the basis functions. 



3.3 Coordinate Transformations 

Let us consider a general coordinate transformation of the 
form x — > x' = (1 + ^)x + e, where $ is a 2 x 2 matrix, 
e = (61,62) is a small displacement. Such a transformation 
can arise for instance from a translation, rotation or from the 
action of gravitational lensing. We assume that the trans- 
formation matrix \E> and the displacement e are small and 
constant across the object. We parametrise the matrix \l/ 
following the gravitational lensing conventions as 



K + 71 
72 + P 



72 - P 

K — 71 



(29) 



where p describes rotations and the convergence k describes 
overall dilatations and contractions. The shear 71 (72) de- 
scribes stretches and compressions along (at 45° from) the 
x-axis. The displacements ei and £2 correspond to transla- 
tions along the x and y-axis, respectively. 

Under this transformation, the intensity /(x) of an ob- 
ject becomes 



/'(x') = /(x(x')) ~ /(x' - *x' - e) 



(30) 



Since we are now considering infinitesimal transformations, 
we can Taylor expand this expression and only keep the 
terms which are first order in \&. After using Equations (jll]) 
and ([|), we find 



/' ~ (1 + pR + nK + jjSj + eifi)f, 



(31) 



where R, K, Si and Ti are the operators generating rota- 
tion, convergence, shears and translations, respectively, and 
where we have used the Einstein summation convention. The 



generators 


are given by 






R = - 


i {Xlf)2 — X2pl) 


= a\d\ — 


a\a,2 


k = - 


i {xipi + X2P2) 


-4' 


~t2 . ~t2 

a[ + a 2 


Si = - 


i {xtpi - X2P2) 


1 (~i2 

2 v 1 


.+2 
— 0, 2 ~ a 


s 2 = - 


i (xip2 + x 2 pi) 


= a\&l - 


a\a,2 


Tj = - 


1 ,. t 


3 = 


~- 1,2. 



-2 ~2\ 

ai - a 2 



(32) 



The rotation generator R is thus simply equal to the angular 
momentum operator in 2-dimensions 



X1P2 — x 2 pi 



(33) 
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up to a factor of —i. Similarly, the translation generator t% 
is simply equal to the the linear momentum operator pi, up 
to the same factor. 

These expressions along with Equation (Q) make it 
easy to compute the effect of these transformations on the 
basis functions B n . For instance, the generator of transla- 
tions along the a;-axis has a matrix representation Ti mn = 
(m|Ti|n) given by 



v2 



V fii + 15 



(34) 



and similarly for the other generators. 

The meaning of the generators can be seen by studying 
their action on the ground state. For instance, it is easy to 
see that under the action of a shear 71, the ground state 1 00) 
becomes 



|00>' ~ (1 + ti&)|00) = |00) + ^= [|20> - |20>] . 

The action of the different transformations on the ground 
states can be calculated in the same way and are shown 
in Figure ^| Clearly, their action is as expected from their 
definition. Since the ground state is circularly symmetric, 
the rotation operator vanishes when applied to 1 00) , i.e. 
i?|00) = 0. More instructively, we can consider the effect 
of R on an asymmetric state like |10). It is also shown in 
the bottom row of the figure. As expected the state 1 10) is 
rotated counter-clockwise by the rotation operator. 

Finite transformations can be produced by exponenti- 
ating the generators. For instance, after a finite rotation by 
an angle p the function / becomes 



(35) 



f' = e pH f= E 



(pR) n 



f, 



(36) 



and similarly for the shear and convergence. 



3.4 Effect of Noise 

In this section, we study the uncertainty induced by noise on 
the basis function decomposition, in the case of correlated 
and uncorrelated background noise, and of Poisson noise. 
The observed intensity of an object is 



/'(x) = /(x)+n(x), 



(37) 



where /(x) is the intrinsic intensity of the object and n(x) is 
the noise. The noise is taken to be unbiased, so that (n(x)) = 
and is characterised by its correlation function rj(x., x') = 
(n(x)n(x )). Here the brackets refer to an ensemble average 
over noise realisations. 

The observed basis coefficients are then = (k,/3\f') 
and are unbiased, i.e. 



</k> = h 



(38) 



where /k = (k, (3\f) axe the intrinsic coefficients. It is then 
easy to show that the covariance matrix cov[/ k , //] = ((/ k — 
(/k))(/i — f° r the observed coefficients is given by 



,/9)»j(x,x')fli(x,/9). (39) 



To be more specific, we first consider the case of ho- 
mogeneous background noise, as can be produced by sky or 



|00> 



(l + e s T 2 )|00> (1+7,S 1 )|00> 



(1 + £,T,)|00> (l+/eK)|00> (l+7 2 S a )|00> 

• • # 




Figure 6. Effect of coordinate transformations on the first two 
shapelet states 1 00) and 1 10) . All coordinate transformations are 
considered: translations, rotations, convergence and shear. They 
are parametrised by 6j, p, K and 7$, respectively, which were all 
assigned a value of 0.3 for the purpose of this figure. Their action 
on any state are easily calculated using the raising and lowering 
operators a J and di. Clearly, the different transformation gener- 
ators act as expected from their definition. 



instrumental backgrounds. If the background noise is uncor- 
related, 7?(x) = a^S^ (x), where a n is the rms noise. As a 
result, the covariance matrix reduces to 

cov[fLJ{] = al5i k , (40) 

where we have used the orthonormality of the basis functions 
(Eq. Q). In this case, the covariance matrix is thus diagonal, 
so that each coefficient is statistically independent. More- 
over, the diagonal elements are all equal. Uncorrelated noise 
thus populates each coefficient equally, and is thus "white" 
as in the case of Fourier transforms. 

In the case of spatially correlated but homogeneous 
noise, the noise correlation function is only a function of 
separation and can thus be written as ??(x — x'). As a result, 
Equation ( pi| ) reduces to the integral of a convolution and 
can thus be written symbolically as 



cov[/ k) /i] = (k,/3|r ? *(l )/ S)>, 



(41) 



in the notation of Equation Q45|) . A convenient way to eval- 
uate this is to decompose rj(x) itself into basis functions and 
then to use the results of §^ below. Spatial correlations in 
the noise thus produces correlations in the coefficients. 

Another case of practical interest is that in which the 
noise is dominated by Poisson shot noise. If the intensities 
are measured in units of photon counts, the noise correlation 



function is then t;(x,x') 
covariance matrix is 



/(x)5 (2) (x-x'). As a result, the 



4/k,/i'] = 



f B {3) 



(42) 
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where (ft, f3, j3) is the 3-product integral defined in 

Equation (^) below, and which is evaluated analytically 
in Paper II. In this case again, the covariance coefficient is 
made non-diagonal by the noise correlation, but is easily 
calculable analytically. 



4 CONVOLUTION 

We now show how shapelets behave under convolutions, an 
operation which often occurs in practice (eg. under the ac- 
tion of PSF, seeing, smoothing, etc) . We start by considering 
convolution by a general kernel in 1-dimensions, and then 
study the special case of smoothing by a gaussian. Finally, 
we treat the 2-dimensional case, and illustrate the results 
with the example of an HST galaxy image. 

4.1 Convolution in 1-Dimension 

Let us first consider the convolution of two arbitrary 1- 
dimensional functions f(x) and g(x). Their convolution h(x) 
can be written as 



h{x) = (J*g)(x) 



dx f(x - x')g{x) 



(43) 



Each function can be decomposed into our basis functions 
with scales a, f3 and 7. These scales are chosen to be 
most convenient in each case. The coefficients are then 
/„ ee (n;a\f), g n = (n;P\g), h„ = {n;j\h}. Our aim is 
to find an expression which relates h n to /„ and g n . Since 
convolution is a bi-linear operation, this relation will be of 
the form 
00 

h n = C nm lfmgi, (44) 

m,l=0 

where the convolution tensor can be written symbolically as 

C nm i(r(,a,P) = (n;j\(m;a) * (l;0)) (45) 

and is a function of the scale lengths. Using the properties 
of the basis functions under Fourier transforms (Eq. ||), it 
is easy to show that the convolution tensor is given by 

C nml ( 7 ,a,P) = (2 7 r)i(-l)V l+m+i Bi 3 j i; ( 7 - 1 ,a- 1 ,/3- 1 ),(46) 

where the 3-product integral is B^^Aai, 0,2, 03) is defined as 

dx Bi(x, ai)B m (x, a,2)B n {x, 03). (47) 



? (3) 



As we show in Paper II, this integral can be easily evaluated 
analytically with a recurrence relation. 

4.2 Smoothing in 1-Dimension 

The special case consisting of smoothing by a gaussian is 
useful in practice. In this case, we let 



g(x) = (2 7 r)"5/3- 1 e 



(48) 



which is normalised so that J dx g(x) = 1. We can then 
write the coefficients for the smoothed function h(x) as 



^ Gnm fm , 



where G„ m (7,a,/3) = J^i Cnmt(j, ot, P)gi is the smooth- 
ing matrix. The gaussian g(x) can be thought as a (non- 
normalised) n — shapelet state of amplitude go — (0; /?[<?}, 
so that Gnm = Cnmogo- Using the generating function for 
Hermite polynomials, one can show that, for the natural 
choice of 7 2 = a 2 + p 2 , the smoothing matrix is given by 



Gnm — 2 2 



(m!/n!)* 



P n a r ' 



( m — n \ 1 
2 )• 



(50) 



for m — n > and even (G nm vanishes otherwise), and where 

— 2 . — 2 1 rj — 2 

lo = a + p 

Figure ^ shows how this analytic formula can be used 
to effici ently smooth a 2-dimensional image (see discussion 
in § |4.3| below). An intuitive feeling for the effect of con- 
volution on the shapelet coefficients can be obtained from 
Figure ^, which graphically shows the smoothing matrix 
G nm (a,f3, 7 = (q 2 +/3 2 )3) for different values of the smooth- 
ing scale j3. As expected, the smoothing matrix approaches 
the identity matrix in the limit of vanishing smoothing scale 
(/3 — » 0). On the other hand, for very large smoothing kernels 
(P — > 00) it reduces to a projection of all the input modes 
m onto the n = output mode. For intermediate scales, the 
smoothing matrix takes the form of a band which rotates 
from the vertical to the horizontal as the smoothing scale 
P is increased. Smoothing thus corresponds to a projection 
of the input modes into output modes of smaller order. The 
high-order modes indeed have oscillations on small scales 
and are thus gradually lost when we increase the smoothing 
scale p. 



4.3 Convolution in 2-Dimensions 

Let us now consider the convolution of two 2-dimensional 
functions, such as 



ft(x) = (/ * g)(x) = / d 2 x /(x - x')ff(x') 



(51) 



We again first decompose each function into shapelet co- 
efficients f n ee (n;a\f), g n = (n;/%), and h n = (n; 7 |/i) 
with shapelet scales a, P and 7 respectively, and where 
n = (n\,n2) as before. Because convolution is bilinear we 
can again relate the convolved to the unconvolved coeffi- 
cients by 



E 



(52) 



where C n ,m,i(7, a, P) is the 2-dimensional convolution ten- 
sor. From the separability of the 2-dimensional basis func- 
tions (see Eq. |l8), it is easy to show that this tensor is equal 
to 

C n , m ,i(7,a,/3) = C„ 1 , mi , !l (7,a,/?)C , „ 2 , m2 , !2 (7,Q,/3), (53) 

where the tensors appearing on the right-hand side are the 
1-dimensional convolution tensor defined in Equation (|l6|). 

We can also consider the special case of smooth- 
ing with a 2-dimensional gaussian. In this case, <?(x) = 
2 

(2ivP 2 )~ 1 e 2^ t which is normalised so that J d 2 x <?(x) = 1. 
The smoothed coefficients are then given by 



(49) 



hn — ^ Gn , m fn 



(54) 
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Figure 7. Illustration of smoothing in Shapelet Space. The orig- 
inal galaxy image (61 X 61 pixels) of Figure wl (shown again in the 
upper-left panel) is smoothed with a gaussian kernel with stan- 
dard deviation = 2 pixels (upper-right panel). The resulting 
image smoothed using shapelets (lower-right panel) is almost in- 
distinguishable from that smoothed using direct convolution in 
real-space (lower-left panel). In shapelet space, smoothing is a 
simple matrix multiplication and can be very efficient when the 
coefficient matrix is sparse, as is the case here (see Figure W). 



where G n ,m(7, «, j3) is the 2-dimensional smoothing matrix. 
It is again easy to show that it is equal to 



Gn.m(7,a,/3) = Gn 1 ,m 1 (7,Q,^)G„ 2 ,m 2 (7,a,/3), 



(55) 



i n te rms of the 1-dimensional smoothing matrix defined in 
§4.2. With the natural choice j 2 = a 2 + (3 2 , it can be eval- 
uated using Equation (pSc]), 

Figure (?] shows the how the galaxy image of Figure ^ 
can be smoothed using our shapelet method. The resulting 
image is indistinguishable from that derived using ordinary 
convolution in real-space (also shown). The shapelet method 
is however computationally very efficient when the coeffi- 
cient matrix is sparse as is the case here (see Figure ^) . The 
effect of smoothing on the shapelet coefficients of this galaxy 
can be seen in Figure ^j. For clarity, the smoothing scale was 
enhanced to /3 = 4 pixels. Clearly, convolution amounts to a 



projection onto the lower order states, as discussed in §4.2 



5 POLAR SHAPELETS 

The cartesian basis functions discussed above are separable 
in the cartesian coordinates x\ and X2- It is also useful to 
construct basis functions which are separable in the polar 
coordinates x and <p. These are eigenstates of the Hamil- 
tonian H and of the angular momentum L simultaneously, 



Figure 8. Graphical representation of the smoothing matrix 
|Gnm| for different size /3 of the smoothing kernel (in units of 
the input shapelet scale a). The horizontal axis corresponds to 
the input (unsmoothed) mode m, while the vertical axis shows the 
output (smoothed) mode n. For small smoothing scales (J3 — » 0) 
the smoothing matrix approaches the identity matrix (upper-left 
panel). For large smoothing scales (/3 — > oo), it approaches a 
projection onto the n = mode (lower-right panel). For inter- 
mediate values, it corresponds to a projection onto lower order 
modes (upper-right and lower- left panels). 



and thus have a number of convenient features. In this sec- 
tion, we show how they can be constructed and study some 
of their properties. 

5.1 Raising and Lowering Operators 

To construct the polar basis functions, we first define the left 
and right lowering operators as (see eg. Cohen- Tannoudji et 
al. 1977) 



1 ,„ , - 1 ,„ 

ai — —= (a i + ia-2) , a r = —= (ai — mi) . 

V 2 v 2 



(56) 



The associated raising operators are a\ and at, respectively. 
The only non-vanishing commutators between these opera- 
tors are 

[ai,al] = [&r, at] = 1. (57) 

The Hamiltonian (Eq. jE)) and angular momentum 
(Eq. [^i)) operators for the 2-dimensional QHO can be then 
be written as 



H = N r + Ni 



L = N r -Ni, 



(58) 



where the left-handed and right-handed number operators 
are naturally defined as 
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Table 1. First few polar Hermite polynomials 
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Figure 9. Shapelet coefficients of the same galaxy (Figure |3|) af- 
ter smoothing with a gaussian kernel. For clarity, the standard 
deviation of the kernel was increased to fi = 4 pixels. By com- 
paring this distribution with that before smoothing (Figure H), it 
is easy to see how convolution amounts to a projection onto the 
lower order shapelet states. 



Ni = ajai, N r = a\.a r 



(59) 



The operators a], a\, ai, and a r can thus be thought of as 
creating and destroying left- and right-handed quanta. 

5.2 Angular Momentum States 

Since the operators iV r and Ni form a complete set of 
commuting observables, their eigenstates \ni,n r ) provide a 
complete basis for our function space. These states are de- 
fined Ni\ni,n r ) = m\ni,n r ) , and similarly for N r , for n; n r 
non-negative integers. They can be constructed by apply- 
ing the raising operators several times on the ground state 
|ni = 0,n 2 = 0} = \ni =0,n r =0) = |0,0>, as 



\n r , ni) 



(at)™- (a?)"' 



In,! 



0,0). 



(60) 



From Equation (^8|), it is easy to see that 

H\m,n r ) = (n r +ni)\n t ,n r ), L\ni,n r ) = (rir— ni)\m, n r ).(61) 

We can therefore relabel these states in terms of their energy 
and angular momentum quantum numbers, n = n r + ni and 
m = n r — m, as 

\n, m) = \m = i(n - m), n r — i(n + m)). (62) 

The angular momentum quantum number takes on the n + 1 
values given by m = — n, —n + 2, . . . , n — 2, n. 
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5.3 Basis Functions 



Using the i-representation of a\ and fit, one can show that 
the basis functions Xni,n r (x,ip) = {x\ni,n r ) for the angular 
momentum states are given by 

Xni ,n r (x,v) = [nm\n r \]-i H n[ , nr (x)e- x2/2 e^-"^, (63) 

where H nk ^ nr (x) are polynomials, which we call 'polar Her- 
mite polynomials'. They can be computed by noting that 
Ho,o{x) — 1 and by using the recursion relation 

l —H k ,i{x) = lH k> i-i-kH k -i,i. (64) 
x 

The diagonal polynomials can be computed using 

Hkk = Hk+i,k-i — x 1 Hk,k-i- (65) 

The first few polar Hermite polynomials are listed in Table |l[ 
They have a number of useful properties. In particular, they 
are symmetric, i.e. Hk,i = Hi,k and their derivative obey 

H' k ,i{x) = kH k -x,i{x)+lH k ,i-x{x) 

= 2xH k ,i(x) - Hk+i,i - Hk,i+i- (66) 

Dimensional polar basis functions can be constructed 

as 

A nu n r (x, ip;/3) = fi~ 1 Xn, l ,n r {P~ 1 x,ip). (67) 
It is easy to check that these are orthonormal, i.e. that 

2T! 



dtp I dx xA ni: „ r (x,<p; f3)A n ^ n , (x,tp; (3) 
..Jo 1 
i i 



{ni,n r ; [3\ni,n r ; [3) = 8, 



(68) 



The radial dependence |X"i,»ir( :r )l °f the first few polar 
shapelet functions is shown in Figure 



5.4 Relation to Cartesian States 

It is useful to relate the angular momentum states \n, m) to 
the cartesian states \n\,ni). Using Equations ( [H^ ) and ( ]60[ ) 
along with the binomial expansion, one can show that the 
transformation matrix between these two bases is given by 



(ni,n 2 |nz,n r ) 



2"3 



i(n r +n;)^n,.-n, 



ni!n 2 ! 



EE- 

ni. =0 n'=0 



. n r \ni \ 
it r 



8ni +!i2,R r +/l; 



This shows that only states with ni + ri2 = n r + rii are mixed. 
The first few |n, m) states are given in terms of \n\n2) states 
in Table § 
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Figure 10. Radial dependence of the first few polar basis func- 
tions \Xn,,n r (x)\. 



Table 2. Angular momentum states \n,m) in terms of the 
cartesian states \n\,n2)- 



: 0,m : 



0} = |0,0> 

> = 
-1) 



l,m=l) = ^[|l,0>+i|0,l>] 
^=[|l,0>-i|0, 1>] 



n - 
n ■■ 

n = 1, m - 

n = 2,m = 2) = \ [|2, 0) + i-y^l, 1) - |0,2>] 
n = 2,m = 0} = ^- [|2,0> + |0,2>] 
n = 2,m = -2) = \ [|2,0> - n/2|l, 1) - |0,2)1 



5.5 Properties 

Because the polar shapelet states are eigenstates of the an- 
gular momentum, they have simple rotational properties. In- 
deed, under a finite rotation by an angle p, the polar states 
transform as 



\n r , 71;)' 



-ipL\ 



,ni), 



(70) 



where we have used the exponentiation (Eq. |56|) of the 
rotation generator R — —iL (Eq. |H) to operate a finite 
rotation. In this basis, finite rotations thus corresponds only 
to a phase factor. 

It is therefore a simple matter to rotate an arbitrary 
function /(x). First, we decompose it into polar shapelet co- 
efficients fn r ,ni — {n r ,ni; /3\f), with an appropriate shapelet 
scale j3. The coefficients = {n r ,ni; /3\f) of the rotated 

function /'(x) are then given simply by 



f 



-i{n r —ni)p 



fn 



(71) 



By contrast, operating a finite rotation in the cartesian 
basis requires an infinite number of applications of the R 
operator (see Eq. |36j) and is thus impractical. On the other 
hand, convolutions do not have simple analytical expressions 
in the polar basis, as they do in the cartesian basis (see §[l| 
and Paper II). The results of §^^, can thus be conveniently 
used to convert from one basis to the other, depending on 
the operation to be performed. 



6 APPLICATIONS 

Now that we have developed the main formalism for 
shapelets, we illustrate the method using images from the 
Hubble Space Telescope (HST). We also discuss several di- 
rect applications of shapelets. 



6.1 Example with HST images 

As an example, we apply the shapelet decomposition method 
to images of galaxies found in the Hubble Deep Field (HDF; 
Williams et al. 1996) , the deepest images observed with the 
Hubble Space telescope (HST). Figure |^ shows the origi- 
nal 61 x 61 image /(x) of one such galaxy (upper-left hand 
panel). Using Equation (^), we first compute the shapelet 
coefficients / n of the galaxy with a shapelet scale f3 = 4 pix- 
els. We then reconstruct the image using Equation ( |23[ ) in- 
cluding coefficients up to a maximum order m +ri2 < n ma x- 
The resulting reconstructed images are shown on the fig- 
ure for different values of n max . As n max is increased, more 
small scale and large scale features are recovered, a s ex - 
pected from the properties of our basis functions (see §2.4). 
For n ma x = 20, the reconstructed image is almost indistin- 
guishable from the original 



Figure |i| gives a graphical representation of the shapelet 
coefficient matrix f ni ,n 2 f° r this image. This can be thought 
as the representation of the galaxy in 'shapelet space'. As 
is apparent on the figure, the coefficient matrix is rather 
sparse. The fact that the odd coefficients are small results 
from the fact that the shapelet center was chosen to be close 
to the centroid of the galaxy. The coefficients are negligibly 
small beyond n\ + ni <; 15, thus explaining why we obtain 
a virtually full reconstruction with n max = 20. 

Because of the sparse nature of the coefficient matrix, 
we can hope to recover the image from only the first few 
largest coefficients. Figure ^| shows the reconstruction of the 
galaxy of Figure | and of two other HDF galaxies, by keeping 
only the A^of coefficients with largest absolute value |/ n |. As 
can be seen from the top two panels, the galaxy image can 
be faithfully recovered with the top AT co f = 60 coefficients, 
yielding a compression factor N p i x /N co { of 62 compared to 
the original image which contained A^pix = 61 x 61 = 3721 
pixels. (The bookeeping required to keep track of selected 
coefficients only requires 1 bit per coefficient, and thus re- 
sults only in a small overhead relatively to this compres- 
sion factor). For the other two galaxies shown in the middle 
and bottom rows, compression factors of about 40-90 are 
achieved. Note that the galaxy in the bottom row has a 
simpler structure and thus affords more compression. 

6.2 Catalogue Archival 

We have seen above that the first few shapelet coefficients 
capture most of the structure of galaxy images and thus al- 
low strong data compressions. This can be very useful for 
upcoming and future large galaxy surveys such as the Sloan 
Digital Sky Survey (SDSS; Gunn & Weinberg 1995), or that 
derived from the planned SNAP mission (Perlmutter et al. 
2001). One can imagine storing the first few shapelet coeffi- 
cients in the catalog, thus both saving storage and conveying 
compactly the shape information of each each galaxy. The 
flux, centroid, major and minor axes and position angle of 
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each galaxy could then be compute d fro m their shapelet co- 
efficients directly (as described in §3.2), thus avoiding the 
need to consider several definitions of magnitudes. Since 
galaxy shapes in different wavelengths are strongly corre- 
lated, the treatment of multi-colour data could be done effi- 
ciently by decomposing differences of the images in different 
pass-bands and again keeping the largest coefficients. The 
resulting catalogue could then also be useful for study of 
galaxy morphology and classification in shapelet space. 



6.3 Modeling the Point-Spread Function 

Several astronomical techniques (eg. high- precision astrom- 
etry and photometry, microlensing, weak lensing, supernova 
searches, etc) require correction for the Point-Spread Func- 
tion (PSF) of the telescope across an image. For instance, 
Alard & Lupton (f998) have developed a technique to take 
the difference of two images convolved with a spatially vary- 
ing PSF. Shapelets provide a convenient correction method 
for the PSF: the PSF shape can be measured at different po- 
sitions in the field using bright stars and then decomposed 
into shapelet coefficients; a 2-dimensional polynomial fit for 
each shapelet coefficient as a function of position can then be 
performed to derive a model of the PSF shape at any point 
(Cf. Alard & Lupton 1998 and Kaiser 2000 who used this ap- 
proach with other sets of basis functions). The convolution 
matrix can then be inverted to compute the shapelet coef- 
ficients of objects prior to convolution. As discussed in 
convolution amounts to a projection onto lower shapelet or- 
der. As a result only low order coefficients can be recovered. 
Another approach consist of fitting the deconvolved shapelet 
coefficients convolved to the PSF model to the data (Cf. 
Kuijken 2000). The analytical properties of shapelets under 
convolution (see §^ and Paper II) greatly facilitate and clar- 
ify the procedure. A detailed study of deconvolution using 
shapelets will be presented in Paper II. 



6.4 Gravitational Lensing 

Gravitational Lensing is a powerful method to directly probe 
the mass of astrophysical objects. In particular, the weak 
coherent distortions induced by lensing on the images of 
background galaxies provide a direct measure of the distri- 
bution of mass in the Universe. This weak lensing method 
is now routinely used to study galaxy clusters, and has re- 
cently been detected in the field (see reviews by Mellier 1999; 
Bartelmann & Schneider; Mellier et al. 2000). Because the 
lensing effect is only of a few percent on large scales, a pre- 
cise method for measuring the shear is required. The original 
methods of Bonnet & Mellier (1995) and of Kaiser, Squires 
& Broadhurst (KSB; 1995) are not sufficiently accurate and 
stable for the upcoming weak lensing surveys. Thus, several 
new methods have been proposed (Kuijken 1999; Rhodes, 
Refregier & Groth 2000; Kaiser 2000). As we briefly de- 
scribe below, the remarkable properties of our basis func- 
tions, make shapelets particularly well suited for providing 
the basis of a new method. 

Let us consider a gal axy with an unlensed intensity 



After decomposing these intensities into our basis functions 
B„(x, (3) (Eq. J22J), this becomes a relation between the 
lensed and the unlensed coefficients given by 



■n rn + 7iSimn)/n 



(73) 



where Si mn = (m|5i|n) is the shear generator matrix given 
in Equation fe2|). The goal for weak lensing is to estimate 
the shear from the shapes of an ensemble of galaxies which 
are assumed to be randomly oriented prior to lensing. In 
the widely used KSB method, this is done by considering 
the effect of lensing on the Gaussian-weighted quadrupole 
moments of the galaxy images. These are exactly equal to 
the ni + 712 = 2 coefficients in our shapelet decomposition. 
In as sense, our method thus generalises this approach and 
capture all the available shape information of galaxies. Be- 
cause the shear matrix is simple in our cartesian basis, we 
can then construct an estimator for the shear by compar- 
ing the distribution of the lensed shapelet coefficients f„ to 
that of a training set /„ for which lensing is known to be 
negligible. This can be done either by constructing a linear 
shear estimator from the observed coefficients or by using a 
Maximum Likelihood technique. In Paper II, we follow the 
first approach and show that shapelets can be used to derive 
precise shear recovery in realistic simulations of deep optical 
images (see Bacon et al. 2001). 

As Luppino & Kaiser (1997) discussed, the shear acts, 
in practice, after the smearing produced by the PSF. To ac- 
count for this, w e ca n first model the PSF across the field, 
as described in § |6.3[ Then one of the methods mentioned 
in this section can be used to derive the deconvolved coeffi- 
cients from the observed convolved ones. Again, a detailed 
study of this deconvolution method and its impact in weak 
lensing measurement will be presented in Paper II. 

Shapelets can also be used for strong lensing applica- 
tions, such as the modeling of cluster or g alax y potentials 
using multiple images and giant arcs. In §3.3, 



we concen- 



/(x). We have shown in §3.3 that under the action of a 



weak shear ji, the lensed intensity is 

f~(i + 1 A)f. 



(72) 



trated mainly on first order distortions, but also mentioned 
that distortions of arbitrary amplitudes can be derived by 
exponentiating the shear and convergence operators (see 
Eq. |$6j). An equivalent method to compute the effect of 
large distortions on the shapelet coefficients is to use the 
analytical expression for rescaling of Appendix |a|. One can 
then model the shape of the lensed object using shapelets 
and explore a large class of lens models efficiently by com- 
puting the lensed image coefficients analytically. Another 
possibility is to also model the gravitational potential of the 
lens using shapelets. 



6.5 Deprojection 

Another important problem in astronomy is that of depro- 
jection. For instance, the 2-dimensional images of galaxies 
and clusters of galaxies observed on the sky are projections 
of their 3-dimensional distributions. One can hope to re- 
construct the 3-dimensional distribution of these systems 
by combining observations at different wavelengths. These 
indeed probe different physical processes, and therefore cor- 
respond to different weighting along the line of sight. Here, 
we show how shapelets can be used to solve this problem. 

To so so, we consider the simple, yet practical ex- 
ample of a cluster of galaxies observed both through its 
X-ray emission (see eg. Sarazin 1988 for a review) and 
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Sunyaev-Zel'dovich (SZ) effect (Sunyaev & Zel'dovich 1972; 
see Birkinshaw 1999, for a review). Cluster deprojection is 
a long standing problem in astrophysics, and has been stud- 
ied by several groups (see the recent solution by Dore et al. 
2001, and references therein). Here, we assume, for simplic- 
ity, that the cluster gas is isothermal. The X-ray emissivity 
of the cluster can then be written as (eg. Sarazin 1988) 



(74) 



X{x,y) ~X I dz p (x,y,z) 



where p is the 3-dimensional electron density of the gas, z 
is the line-of-sight coordinate, and Xq is a constant which 
depends on the wavelength of observation, the gas temper- 
ature, and the distance to the cluster. The comptonisation 
parameter Y(x, y) from the SZ effect can be observed as tem- 
perature anisotropies of the Cosmic Microwave Background 
(see review by Birkinshaw 1999). It is proportional to the 
electron pressure integrated along the line of sight, and thus 
be written, for an isothermal cluster as 



Y(x,y) ~Y dz p(z,y,z) 



(75) 



where Yo is a constant which again depends on the wave- 
length of observation, the gas temperature, and distance to 
the cluster. Our goal is to reconstruct the three-dimensional 
gas density p(x,y,z) from measurements of X(x,y) and 
Y(x,y). 

For this purpose, let us decompose these two observed 
images into 2-dimensional shapelets 

as X(x,y) = J2 nin2 X "i"2 B n 1 n 2 (x,y), and Y(x,y) = 
X/n n ^ni"2^nin2 (x , y) ■ We choose the same shapelet scale 
f3 for X, Y and p and thus drop it to sim plify the no- 
tation. In analogy with the discussion in 
also define 3-dimensional basis functions B 
B ni (x)B n2 (y)B„ 3 (z) as products of three 1-dimensional 
shapelets. This allows us to decompose the 3-dimensional 
gas density distribution as 



3.1 



^ni no n r - 



we can 
(x,y,z) = 



p{x,y,z) 



E 

rc.i ,n 2 ,n 3 



/?n i n 2 n 3 B n 1 n 2 n 3 ?/, z). 



(76) 



Using the properties of the shapelet basis functions, it is 
then easy to show that the shapelet coefficients for the X- 
ray emissivity can be written as 



X n 



= x 



E 



B 



(3) 



(3) 



' PmPn 



(77) 



(3) 

where B„', is the ubiquitous 3-product integral defined in 
Equation (H7J), and m = (mi, mi, W.3) in this context. Sim- 
ilarly, the coefficients for the comptonisation parameter can 
be written as 



™3 



l|n 3 )p n , 



(78) 



where (l|n) is the integral defined in Equation (^). The di- 
rect approach, which consists in solving these two equations 
of the desired coefficients p n , is probably difficult in practice. 
A more convenient approach is to derive an estimate for p n 
by x 2 -htting these coefficients to the observables X ni „ 2 and 
Y ni n 2 taking into account the noise in each measurement. 
The \ 2 procedure also produces the covariance matrix for 
the coefficients p n , and thus allows us to study any degen- 
eracy present in the deprojection. This is greatly facilitated 



in practice by the analytic forms for (l|n) (Eq. E7J) and for 
Bnmi ( see P a P er II) , and the fact that the fitted model is 
linear in its output parameters p n . 

Note that our method is fully general and does not as- 
sume that the cluster distribution has any specific form. In 
particular, it could be particularly useful if the SZ observa- 
tions are performed with an interferometer as is the case for 
recent measurements (see Carlstrom et al. 1999, and refer- 
ence therein). In this case, the interferometer yields a mea- 
surement of the Fourier transform of Y(x,y), and can thus 
make use of the dual properties of our shapelet functions 
under Fourier transforms (Eq. Q). A more thorough study 
of the deprojection using shapelets is left to future work. 
A study of the use of shapelets for reconstructing images 
with interferometers will be presented in Chang & Refregier 
(2001). 



7 CONCLUSIONS 

We have described and developed a new method for 
analysing images. It is based on the decomposition of the ob- 
jects in the image into a series of basis functions of different 
shapes, or 'shapelets'. The method is fully linear and uses a 
number of powerful properties of the basis functions. In par- 
ticular, we showed that Hermite basis functions have simple 
analytic properties under convolution, noise, rotations, dis- 
tortions, and rescaling. These functions are eigenfunctions 
of the QHO and thus allow us to use the formalism devel- 
oped for this problem. For instance, we showed that trans- 
formations such as translations, rotations, shears and dilata- 
tions can be expressed as simple combinations of the rais- 
ing and lowering operators. Another remarkable property of 
these functions is that they are (up to a rescaling) their own 
Fourier transforms. This is a unique property, which stems 
from the special symmetry of the QHO Hamiltonian. We de- 
rived analytical expressions for the flux, centroid and radius 
of the object, from its shapelet coefficients. We also con- 
structed polar shapelets which give the explicit rotational 
properties of the object. 

It is interesting to compare our method to the wavelet 
method (see review by Stark, Murtagh & Bijaoui 1998). In 
this latter method, the image is decomposed into a sum of 
basis functions located on a grid across the image. The ba- 
sis functions are taken to have a range of sizes, but have all 
the same shape. Wavelets are thus ideal to decompose an 
image into different scales, which can then be analysed sep- 
arately. Our method, on the other hand models the image as 
a collection of discrete objects of arbitrary shapes and sizes. 
It is therefore particularly well adapted to the treatment 
of astronomical images, which are typically composed of a 
superposition of compact disjoint objects. The two meth- 
ods can thus be thought as complementary. For instance, 
one can use wavelets to remove large-scale background vari- 
ations, and to search and detect objects in the image. The 
resulting object catalog can then be used as the input to the 
shapelet method, which will then characterise the shape of 
each object in detail. 

Our method potentially has a wide range of applica- 
tions. It can be viewed as a new representation of images 
which makes object shapes easy to study and modify. For 
instance, we applied our method to galaxy images found in 
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the HDF and showed how they could be well represented 
with a small number of shapelet coefficients. This can be 
used to compress galaxy images by a factor of 40-90, and 
could thus have important applications for galaxy archival. 
We also discussed several direct applications of shapelets to 
measurements of gravitational lensing, and the problems of 
de-projection and PSF correction. Other applications to be 
explored are that of multi-colour shapelets and of the study 
of galaxy morphology and classification using shapelets. Our 
original motivation for developing this method was to find a 
robust and precise method to measure weak lensing distor- 
tions in the presence of a PSF. The application of shapelets 
to this problem and to the general problem of deconvolu- 
tion will be presented in detail in Paper II. The application 
of shapelets to image reconstructions from interferometric 
observations will be presented in Chang & Refregier 2001. 
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APPENDIX A: RESCALING 

In this appendix, we show how we can easily operate a 
change of scale for the decomposition of a function in 
1-dimension. This is convenient for finding the optimal scale 
for a given function. In addition, such a change of scale 
occurs when a 2-dimensional image is distorted by gravita- 
tional lensing. 

Let us consider a function f(x) = (x\f) which we de- 
compose (as in Eq. [9) into two sets of basis functions with 
scales fii and 02 as 



f(x) = ^(n;/3i|/) B n (x;0i) = ^{n;0 2 \f} B n (x; /3 2 )(Al) 
The coefficients in each basis are related by 



\ " 



(m;/3i|/>= > (ni;/3i\n 2 ;02){n2;02\f). 



n 2 =0 



(A2) 



Using the generating function of Hermite polynomials, one 
can show that the transformation matrix is given by 



min(roi ,ri2 ) 

<n i; /?i|n 2 ;/3 2 ) = ^ ^ 
(=0 



(m!n 2 !) - 



where 

, _ Pi - Pi , 



;n(m,n 2 ,0 (y) 
2/3 1 j3 2 



b 2 + \ (A3) 



(A4) 



and Il(ni,n 2 ,Z) is equal to 1 if m, n 2 and I are all odd 
or all even, and is equal to otherwise. In the limiting 
case where 0i =02, the transformation matrix reduces to 
(n\\ 0\ |n 2 ; 2 ) = S ni>rl2 , in agreement with the orthonormal 
properties of the basis functions (Eq. [^). 
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